function [AB, AL, AC, AR, AT, S, dx, dy, dv, dc, siga, sct] = ...
    twoDcoefMG_corr( dat, abdat, numg, numm, xcm, xfm, ycm, yfm, mt, src, it, BC, ...
    IBSL, IBSR, IBSB, IBST)

% function A = twoDcoef( dat, numg, numm, xcm, xfm, ycm, yfm, mt, src, it, BC,
%    IBSL, IBSR, IBSB, IBST)
%   This function produces the coefficient matrix and source for solution
%   of the two-dimensional diffusion equation.

% expand the bc's
BCL = BC(1);
BCR = BC(2);
BCB = BC(3);
BCT = BC(4);

% define data components
dc = zeros(numm,numg); sa=dc; ab = dc;
sc = zeros(numm,numg,numg);
for m = 1:numm
    for g = 1:numg %g+(m-1)*g
        dc(m,g) = dat(numg*(m-1)+g,1);   % diffusion coefficient
        sa(m,g) = dat(numg*(m-1)+g,2);   % removal cross-section
        ab(m,g) = abdat(numg*(m-1)+g,1); % absorption cross-section
        ns(m,g) = dat(numg*(m-1)+g,3); % fission cross-section
        for gg = 1:numg
            sc(m,g,:) = dat(g+(m-1)*g,4:end);
        end
    end
end
% number of fine meshes in x and y coordinates
N = sum(xfm);
M = sum(yfm);
% number of coarse meshes in x and y coordinates
CN = length(xfm);
CM = length(yfm);
% compute dx and dy vectors
dx = zeros(N,1);
dy = zeros(M,1);
% coarse mesh index for fine mesh
cix=zeros(N,1);
ciy=zeros(M,1);
j = 0;
for i = 1:CN
    dx( (j+1):(j+xfm(i)) ) = (xcm(i+1)-xcm(i))/xfm(i);
    cix( (j+1):(j+xfm(i)) )=i;
    j = sum(xfm(1:i));
end
j = 0;
for i = 1:CM
    dy( (j+1):(j+yfm(i)) ) = (ycm(i+1)-ycm(i))/yfm(i);
    ciy( (j+1):(j+yfm(i)) )=i;
    j = sum(yfm(1:i));
end
% mesh edge average volume
dv = zeros( (N+1)*(M+1), 1 );

lenk = (N+1)*(M+1);
AC   = zeros( lenk, numg );
AL   = zeros( lenk-1, numg );
AR   = AL;
AB   = zeros( lenk-N-1, numg );
AT   = AB;
S    = zeros( (N+1)*(M+1), numg ); siga = S;
sct  = zeros( (N+1)*(M+1), numg, numg );

for g = 1:numg % group loop
    
    % interior coefficients
    for i = 2:N
        for j = 2:M

            m11 = mt( cix(i-1), ciy(j-1) );  % 12  |  22
            m21 = mt( cix(i  ), ciy(j-1) );  % ____|____ (mat location)
            m12 = mt( cix(i-1), ciy(j  ) );  %     |
            m22 = mt( cix(i  ), ciy(j  ) );  % 11  |  21
            
            k = i+(j-1)*( N+1 ); % 10 being the number of x ints

            alft = -0.25*(dc(m11,g)+dc(m12,g))*(dy(j-1)+dy(j))/dx(i-1);
            argt = -0.25*(dc(m21,g)+dc(m22,g))*(dy(j-1)+dy(j))/dx(i);
            abot = -0.25*(dc(m11,g)+dc(m21,g))*(dx(i-1)+dx(i))/dy(j-1);
            atop = -0.25*(dc(m12,g)+dc(m22,g))*(dx(i-1)+dx(i))/dy(j);
            a    = 0.25*(dx(i-1)*dy(j-1)*sa(m11,g) + ...
                         dx(i)  *dy(j-1)*sa(m21,g) + ...
                         dx(i-1)*dy(j)  *sa(m12,g) + ...
                         dx(i)  *dy(j)  *sa(m22,g) ) - ...
                         (alft+argt+abot+atop);       
                     
            AL( k-1, g )   = alft;
            AR( k, g )     = argt;
            AB( k-N-1, g ) = abot;
            AT( k, g )     = atop;
            AC( k, g )     = a;  % ac ad  bc bd
            
            sct(k,g,:) = 0.25*( dx(i-1)*dy(j-1)*sc(m11,g,:) + ...
                                dx(i)  *dy(j-1)*sc(m21,g,:) + ...
                                dx(i-1)*dy(j)  *sc(m12,g,:) + ...
                                dx(i)  *dy(j)  *sc(m22,g,:) );

            S( k, g )  = 0.25*( dx(i-1)*dy(j-1)*ns(m11,g) + ...
                                dx(i)  *dy(j-1)*ns(m21,g) + ...
                                dx(i-1)*dy(j)  *ns(m12,g) + ...
                                dx(i)  *dy(j)  *ns(m22,g) );
                            
            siga(k,g)  = 0.25*( dx(i-1)*dy(j-1)*ab(m11,g) + ...
                                dx(i)  *dy(j-1)*ab(m21,g) + ...
                                dx(i-1)*dy(j)  *ab(m12,g) + ...
                                dx(i)  *dy(j)  *ab(m22,g) );
                            
            if g == 1
             dv( k, 1 ) = 0.25*( dx(i-1)*dy(j-1) + ...
                                 dx(i)  *dy(j-1) + ...
                                 dx(i-1)*dy(j)   + ...
                                 dx(i)  *dy(j) );
            end
        end
    end
    
    % left and right edges
    for j = 2:M

        i = 1; % LEFT EDGE
        k = i+(j-1)*( N+1 );

        if BCL==1
            beta=1;
        else
            beta=0;
        end
         
        m21 = mt( cix(i  ), ciy(j-1) );  % ____|____ (mat location)
        m22 = mt( cix(i  ), ciy(j  ) );  % 11  |  21

        argt = -0.25*(dc(m21,g)+dc(m22,g))*(dy(j-1)+dy(j))/dx(i);
        abot = -dc(m21,g)*dx(i)/dy(j-1);
        atop = -dc(m22,g)*dx(i)/dy(j);

        a    = 0.25*(dx(i)  *dy(j-1)*sa(m21,g) + ...
                     dx(i)  *dy(j)  *sa(m22,g) ) - ...
                     (argt+abot+atop) + ...
               0.5*(1-beta)/(1+beta);          

        AR( k, g )     = argt;
        AB( k-N-1, g ) = abot;
        AT( k, g )     = atop;
        AC( k, g )     = a;  % ac ad  bc bd

        S( k, g )  = 0.25*( dx(i)*dy(j-1)*ns(m21,g) + ...
                            dx(i)*dy(j)  *ns(m22,g) );
        sct(k,g,:) =  0.25*( dx(i)  *dy(j-1)*sc(m21,g,:) + ...
                             dx(i)  *dy(j)  *sc(m22,g,:) );                
        if BCL == 2 % add the incident source
            S(k,g) = S(k,g) + 0.5*(IBSL(ciy(j-1),g)+IBSL(ciy(j),g));
        end
        if g==1                
            dv( k, 1 ) = 0.25*( dx(i)  *dy(j-1) + ...
                                dx(i)  *dy(j) );
        end            
        
        i = N+1; % RIGHT EDGE
        k = i+(j-1)*( N+1 );    
        
        if BCR==1
            beta=1;
        else
            beta=0;
        end
        
        m11 = mt( cix(i-1), ciy(j-1) );  % 12  |  22
        m12 = mt( cix(i-1), ciy(j  ) );  %     |

        alft = -0.25*(dc(m11,g)+dc(m12,g))*(dy(j-1)+dy(j))/dx(i-1);
        abot = -(dc(m11,g))*(dx(i-1))/dy(j-1);
        atop = -(dc(m12,g))*(dx(i-1))/dy(j);
        a    = 0.25*(dx(i-1)*dy(j-1)*sa(m11,g)  + ...
                     dx(i-1)*dy(j)  *sa(m12,g)) - ...
                     (alft+abot+atop) + ...
               0.5*(1-beta)/(1+beta);
           
        AL( k-1, g )   = alft;
        AB( k-N-1, g ) = abot;
        AT( k, g )     = atop;
        AC( k, g )     = a;  
        
        S( k, g )  = 0.25*( dx(i-1)*dy(j-1)*ns(m11,g) + ...
                            dx(i-1)*dy(j)  *ns(m12,g) );
        sct(k,g,:) =  0.25*( dx(i-1)*dy(j-1)*sc(m11,g,:) + ...
                             dx(i-1)*dy(j)  *sc(m12,g,:) );                        
        if BCR == 2 % add the incident source
            S(k,g) = S(k,g) + 0.5*(IBSR(ciy(j-1),g)+IBSR(ciy(j),g));
        end
        if g==1                
            dv( k, 1 ) = 0.25*( dx(i-1)  *dy(j-1) + ...
                                dx(i-1)  *dy(j) );
        end    
        
    end
    
    % top and bottom edges 
    for i = 2:N
        
        j = 1; % BOTTOM EDGE
        k = i+(j-1)*( N+1 );
        
        if BCB==1
            beta=1;
        else
            beta=0;
        end
        
        m12 = mt( cix(i-1), ciy(j  ) );  %     |
        m22 = mt( cix(i  ), ciy(j  ) );  % 11  |  21

        alft = -(dc(m12,g))*(dy(j))/dx(i-1);
        argt = -(dc(m22,g))*(dy(j))/dx(i);
        atop = -0.25*(dc(m12,g)+dc(m22,g))*(dx(i-1)+dx(i))/dy(j);
        a    = 0.25*(dx(i-1)*dy(j)  *sa(m12,g) + ...
                     dx(i)  *dy(j)  *sa(m22,g) ) - ...
                     (alft+argt+atop) + ...
                     0.5*(1-beta)/(1+beta);           
        AL( k-1, g )   = alft;
        AR( k, g )     = argt;
        AT( k, g )     = atop;
        AC( k, g )     = a;  % ac ad  bc bd

        sct(k,g,:) =  0.25*( dx(i-1)*dy(j)  *sc(m12,g,:) + ...
                             dx(i)  *dy(j)  *sc(m22,g,:) );
        
        S( k, g )  = 0.25*( dx(i-1)*dy(j)  *ns(m12,g) + ...
                            dx(i)  *dy(j)  *ns(m22,g) );
                   
        if BCB == 2 % add the incident source
            S(k,g) = S(k,g) + 0.5*(IBSB(cix(i-1),g)+IBSB(cix(i),g));
        end
        if g==1                
            dv( k, 1 ) = 0.25*( dx(i)  *dy(j) + ...
                                dx(i-1)*dy(j) );
        end    
 
        j = M+1; % TOP EDGE
        k = i+(j-1)*( N+1 );
        
        if BCT==1
            beta=1;
        else
            beta=0;
        end
        
        m11 = mt( cix(i-1), ciy(j-1) );  % 12  |  22
        m21 = mt( cix(i  ), ciy(j-1) );  % ____|____ (mat location)

        alft = -(dc(m11,g))*(dy(j-1))/dx(i-1);
        argt = -(dc(m21,g))*(dy(j-1))/dx(i);
        abot = -0.25*(dc(m11,g)+dc(m21,g))*(dx(i-1)+dx(i))/dy(j-1);
        a    = 0.25*(dx(i-1)*dy(j-1)*sa(m11,g) + ...
                     dx(i)  *dy(j-1)*sa(m21,g) ) - ...
                     (alft+argt+abot) + ...
                     0.5*(1-beta)/(1+beta);     
        AL( k-1, g )   = alft;
        AR( k, g )     = argt;
        AB( k-N-1, g ) = abot;
        AC( k, g )     = a;  % ac ad  bc bd

        sct(k,g,:) =  0.25*( dx(i-1)*dy(j-1)*sc(m11,g,:) + ...
                             dx(i)  *dy(j-1)*sc(m21,g,:) );
                         
        S( k, g )  = 0.25*( dx(i-1)*dy(j-1)*ns(m11,g) + ...
                            dx(i)  *dy(j-1)*ns(m21,g) );
                   
        if BCT == 2 % add the incident source
            S(k,g) = S(k,g) + 0.5*(IBST(cix(i-1),g)+IBST(cix(i),g));
        end
        if g==1                
            dv( k, 1 ) = 0.25*( dx(i-1)*dy(j-1) + ...
                                dx(i)  *dy(j-1) );
        end    

    end
    
    % finally, the corners

    % BOTTOM LEFT ---------------------------------------------------------
    i=1; j=1;
    
    if BCL==1
        beta1=1;
    else
        beta1=0;
    end
    if BCB==1
        beta2=1;
    else
        beta2=0;
    end
    
    m22 = mt( cix(i  ), ciy(j  ) );  % 11  |  21
    k = i+(j-1)*( N+1 ); % 10 being the number of x ints

    argt = -(dc(m22,g))*(dy(j))/dx(i);
    atop = -(dc(m22,g))*(dx(i))/dy(j);
    a    = 0.25*(dx(i)*dy(j)*sa(m22,g) ) - ...
                 (argt+atop)+...
                 0.5*(1-beta1)/(1+beta1) + 0.5*(1-beta2)/(1+beta2);           
    AR( k, g )     = argt;
    AT( k, g )     = atop;
    AC( k, g )     = a;  % ac ad  bc bd

    sct(k,g,:) =  0.25*( dx(i)*dy(j)*sc(m22,g,:) );
    S( k, g )  =  0.25*( dx(i)*dy(j)*ns(m22,g)   );
    if BCL == 2 % add the (left) incident source
        S(k,g) = S(k,g) + 0.5*(IBSL(ciy(j),g));  % NOTE! Not sure about 1/2*********
    end
    if BCB == 2 % add the (bottom) incident source
        S(k,g) = S(k,g) + 0.5*(IBSB(cix(i),g));
    end    
    if g==1                
        dv( k, 1 ) = 0.25*( dx(i)*dy(j) );
    end 

    % TOP LEFT ------------------------------------------------------------
    i=1; j=M+1; 
    
    if BCL==1
        beta1=1;
    else
        beta1=0;
    end
    if BCT==1
        beta2=1;
    else
        beta2=0;
    end
    
    m21 = mt( cix(i  ), ciy(j-1) );  % 11  |  21
    k = i+(j-1)*( N+1 ); % 10 being the number of x ints

    argt = -(dc(m21,g))*(dy(j-1))/dx(i);
    abot = -(dc(m21,g))*(dx(i))/dy(j-1);
    a    = 0.25*(dx(i)  *dy(j-1)  *sa(m21,g) ) - ...
                 (argt+abot)+...
                 0.5*(1-beta1)/(1+beta1) + 0.5*(1-beta2)/(1+beta2);           
    AR( k, g )     = argt;
    AB( k-N-1, g ) = abot;
    AC( k, g )     = a;  % ac ad  bc bd

    sct(k,g,:) =  0.25*( dx(i)*dy(j-1)*sc(m21,g,:) );
    S( k, g )  =  0.25*( dx(i)*dy(j-1)*ns(m21,g)   );
    if BCL == 2 % add the (left) incident source
        S(k,g) = S(k,g) + 0.5*(IBSL(ciy(j-1),g));  % NOTE! Not sure about 1/2*********
    end
    if BCT == 2 % add the (bottom) incident source
        S(k,g) = S(k,g) + 0.5*(IBST(cix(i),g));
    end    
    if g==1                
        dv( k, 1 ) = 0.25*( dx(i)*dy(j-1) );
    end 
    
    % BOTTOM RIGHT --------------------------------------------------------
    i=N+1; j=1; 
    if BCR==1
        beta1=1;
    else
        beta1=0;
    end
    if BCB==1
        beta2=1;
    else
        beta2=0;
    end
    
    m12 = mt( cix(i-1), ciy(j  ) );  % 11  |  21
    k = i+(j-1)*( N+1 ); % 10 being the number of x ints

    alft = -(dc(m12,g))*(dy(j))/dx(i-1);
    atop = -(dc(m12,g))*(dx(i-1))/dy(j);
    a    = 0.25*(dx(i-1)*dy(j)*sa(m12,g) ) - ...
                 (alft+atop)+...
                 0.5*(1-beta1)/(1+beta1) + 0.5*(1-beta2)/(1+beta2);           
    AL( k-1, g )   = alft;
    AT( k, g )     = atop;
    AC( k, g )     = a;  % ac ad  bc bd

    sct(k,g,:) =  0.25*( dx(i-1)*dy(j)*sc(m12,g,:) );
    S( k, g )  =  0.25*( dx(i-1)*dy(j)*ns(m12,g)   );
    if BCL == 2 % add the (left) incident source
        S(k,g) = S(k,g) + 0.5*(IBSR(ciy(j),g));  % NOTE! Not sure about 1/2*********
    end
    if BCB == 2 % add the (bottom) incident source
        S(k,g) = S(k,g) + 0.5*(IBSB(cix(i-1),g));
    end    
    if g==1                
        dv( k, 1 ) = 0.25*( dx(i-1)*dy(j) );
    end 
    
    % TOP RIGHT -----------------------------------------------------------
    i=N+1; j=M+1;
    
    if BCR==1
        beta1=1;
    else
        beta1=0;
    end
    if BCT==1
        beta2=1;
    else
        beta2=0;
    end
    
    m11 = mt( cix(i-1), ciy(j-1) );  % 11  |  21
    k = i+(j-1)*( N+1 ); % 10 being the number of x ints

    alft = -(dc(m11,g))*(dy(j-1))/dx(i-1);
    abot = -(dc(m11,g))*(dx(i-1))/dy(j-1);
    a    = 0.25*(dx(i-1)  *dy(j-1)  *sa(m12,g) ) - ...
                 (alft+atop)+...
                 0.5*(1-beta1)/(1+beta1) + 0.5*(1-beta2)/(1+beta2);           
    AL( k-1, g )   = alft;
    AB( k-N-1, g ) = abot;
    AC( k, g )     = a;  % ac ad  bc bd

    sct(k,g,:) =  0.25*( dx(i-1)*dy(j-1)*sc(m11,g,:) );
    S( k, g )  =  0.25*( dx(i-1)*dy(j-1)*ns(m11,g)   );
    if BCL == 2 % add the (left) incident source
        S(k,g) = S(k,g) + 0.5*(IBSR(ciy(j-1),g));  % NOTE! Not sure about 1/2*********
    end
    if BCB == 2 % add the (bottom) incident source
        S(k,g) = S(k,g) + 0.5*(IBST(cix(i-1),g));
    end    
    if g==1                
        dv( k, 1 ) = 0.25*( dx(i-1)*dy(j-1) );
    end 


end % group loop

